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For the first time, we determine the complete spin- and momentum-dependent 
order v 2 corrections to the static interquark potential from simulations of QCD 
in the valence quark approximation at inverse lattice spacings of 2-3 GeV. 
A new flavor dependent correction to the central potential is found. We re- 
port a r~ 2 contribution to the long range spin-orbit potential V{. The other 
spin-dependent potentials turn out to be short ranged and can be well under- 
stood by means of perturbation theory. The momentum-dependent potentials 
' qualitatively agree with minimal area law expectations. In view of spectrum 



calculations, we discuss the matching of the effective nonrelativistic theory 
to QCD as well as renormalization of lattice results. In a first survey of the 
resulting bottomonia and charmonia spectra we reproduce the experimental 
levels within average errors of 12.5 MeV and 22 MeV, respectively. 
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I. INTRODUCTION 



Quarkonia spectroscopy provides a wealth of information and thus constitutes an im- 
portant observational window to the phenomenology of confining quark interactions. It has 
been known for a long time that purely phenomenological or QCD inspired potential models 
offer a suitable heuristic framework to understand the empirical charmonium (J/i/j) an d 
bottomonium (T) spectra 

On a more fundamental level, one would prefer to start out from the basic QCD La- 
grangian to solve the heavy quarkonia bound state problem. NRQCD || offers a systematic 
way to solve this problem by direct extraction of bound state masses from effective nonrela- 
tivistic lattice Lagrangians, which approximate the QCD Lagrangian to a given order in the 
quark velocity, v. Considerable success has been achieved recently in determining quarkonia 
spectra within this approximation to QCD ||. 

Here, we follow a complementary strategy: instead of separately computing the spectral 
properties of individual mesonic states, we integrate out the gauge background and directly 
determine the underlying quantum mechanical two-particle Hamiltonian. Once QCD bind- 
ing problems are recast into this form, spectra, wave functions and decay constants for 
arbitrary (sufficiently large) quark masses and quantum numbers can easily be obtained. 
Results can either be confronted with experiment or compared to predictions from lattice 
NRQCD. 

In the limit of infinite quark mass, the Born-Oppenheimer approximation is applicable 
and, after integrating out the gauge degrees of freedom, QCD binding problems become 
nonrelativistic. The static interaction potential can be computed directly from the QCD 
Lagrangian on the lattice. Within the present study, we find the average velocity between the 
sources to be (v 2 ) ~ 0.27 and (v 2 ) ~ 0.08 for the charmonium and bottomonium ground- 
states, respectively. This leads us to expect that the phenomenological potentials within 
those models, which have been optimized to reproduce empirical spectra, should deviate 
by substantial 0(v 2 ) corrections from the static potential as predicted by QCD; at realistic 
quark masses such corrections, that are also required to obtain hyperfine splittings, cannot 
be neglected. Therefore, we have to take corrections to the static limit into account. 

The Hamiltonian that we derive is equivalent to the QCD Lagrangian up to 0(v 2 ). It 
includes the spin-dependent (SD) terms derived by Eichten, Feinberg and Gromes J7|||, 
the momentum-dependent (MD) corrections derived by Barchielli, Brambilla, Montaldi and 
Prosperi (BBMP) and one-loop radiative corrections from matching the effective theory 



to the full theory at a scale \i that, in general, differs from the heavy quark mass, m ||10|| . It 
can be parametrized in terms of seven independent scalar functions of the quark separation 
(the potentials). These will be computed nonperturbatively on the lattice. 

The static potential has been determined with high accuracy in the valence quark 
(quenched) approximation to QCD [TTL|r5| and, more recently, in full QCD with two dynam- 
ical flavors of light Wilson sea quarks |14| . First attempts to compute relativistic corrections 
have been made in the mid eighties for SU(2) and SU(3) gauge theory JToLjTS] and have been 
extended to QCD with sea quarks in Refs. [fT9| , p0]| . 

In view of the general interest in the Hamiltonian formulation of the meson binding 
problem, renewed effort should be made to unravel the structure of the SD potentials and 
other 0{v 2 ) corrections. Recently, we presented improved techniques for computation of 
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SD corrections and tested them successfully on SU(2) gauge theory [pl[| . Here, we shall 
apply these methods to the physically relevant SU(3) gauge theory. We will extend the 
SU(2) investigation by inclusion of MD potentials and relativistic corrections to the central 
potential, and subsequently determine quarkonia levels. 

We wish to emphasize that the method presented does not rely on any other approxima- 
tions than truncating the QCD Lagrangian at second order in the quark velocity (apart from 
the valence quark approximation). However, the Schrodinger-Pauli approach to heavy quark 
binding problems suffers from the same difficulties as NRQCD, namely (a) the error involved 
in truncating the expansion at a finite order in v, (b) uncertainties in the matching of the 
effective Hamiltonian to QCD and (c) renormalization of lattice results. While we manage 
to solve the latter problem in a satisfactory way, we have to rely on one-loop perturbation 
theory for the matching of the nonrelativistic Hamiltonian to QCD. Systematic errors from 
the 0(v 2 ) approximation as well as from the uncertainty in the matching constants (that 
can be reduced order by order in perturbation theory) are estimated. 

Since NRQCD to order v 2 (or v 4 , depending on the labeling conventions used) is based on 
the same Lagrangian, it is worthwhile to compare the two approaches. While NRQCD can in 
principle be generalized to any order in v, the Schrodinger-Pauli approach is only valid up to 
order v 2 . Also, we cannot treat heavy-light systems. In NRQCD the zero point energy can be 
fixed by measuring the dispersion relation while in our approach only properties of particles 
at rest can be studied. The clear advantage of the method presented here is that with 
one simulation only, we easily obtain all spectral properties (including arbitrary excitations) 
for any (sufficiently heavy) quark mass. From a two-body Hamiltonian formulation of the 
problem, the effect of individual terms on the spectrum becomes immediately apparent, 
and a transparent understanding of the anatomy of the underlying interaction mechanism is 
obtained. The potentials are protected by the global Z 3 symmetry [32] from finite size effects, 
contrary to NRQCD wave functions and masses, such that we can determine the potentials 
for the r-range required, even for broad excited state wave functions, on relatively small 
spatial lattice volumes. 

The article is organized as follows: in Sec. |T|, we introduce the Hamiltonian and present 
definitions of the potentials that are suitable for lattice evaluation. Moreover, we include 
theoretical expectations on the form of the potentials. Sec. |H] contains simulation de- 
tails and lattice specific techniques wherever they differ from our SU(2) investigation [2T 



The renormalization of lattice operators and the matching procedure between the effective 
nonrelativistic theory and QCD are discussed in Sec. [TV]. The resulting SU(3) potentials 
are presented in Sec. |V[ Promising results on charmonium and bottomonium spectra are 
obtained and discussed in Sec. [VT|, before we conclude. 



II. THE HEAVY QUARK POTENTIAL 
A. Hamiltonian formulation of the meson binding problem 



In Ref. ||21|| , we restricted ourselves to an evaluation of SD corrections to the static 
potential. Since we are going to include the complete 0(v 2 ) corrections into the present 
study and aim to predict quarkonia properties, we find it worthwhile to briefly sketch some 
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details of the derivation of the Hamiltonian. The SD and MD parts as well as relativistic 
corrections to the central potential have been derived during the eighties 0HH • The matching 
problem between QCD and the effective Hamiltonian has been sorted out to one-loop order 
for the SD terms recently |T(| and we extend this to the remaining corrections. 

It is instructive to start at 0(v ), before proceeding to the 0(v 2 ) Hamiltonian. To this 
order, the heavy quark propagator S(x, y) = Q(x)Q*(y) of a quark with mass m obeys the 
evolution equation in an external gauge field AJ^, 



( D 2 \ 

- d 4 S(x,y) = UgA A + m - — I S(x,y), 



(1) 



where _D M denotes the covariant derivative. To O(v ), the solution to the initial value 
problem, 



S(x,y)\ 



£4=2/4 



<ffx 



is given by, 



S(x,y) = U(x;x 4 ,y 4 ) Texp 



J/4 



dt m 



^'4 



2m 



5 3 (x-y), 



(2) 



(3) 



where T denotes the time ordering operator. t/(x; x 4 , y 4 ) is the static propagator of a 
quark, traveling from the point (x, x 4 ) to (x, y 4 ), and consists of the corresponding temporal 
Schwinger line times the factor exp (—E t), with r = y 4 — x 4 . E (fi) represents the static 
quark self-energy that diverges like /z/ln/z with the cut-off scale /1. 

By combining two static propagators into a Wilson loop, one can determine the potential 
Vo(r) between two static sources, separated by a distance r, in the limit of large Euclidian 
times, 



(W(r, r)) oc exp (-V (r)r) , (1 



00 



(4) 



Note that the potential contains the static quark self-energies. In order to obtain the spec- 
trum of mesonic heavy quark bound states, the Schrodinger equation (in the cm. frame), 



Hip n i m (r) 

can be solved, where the Hamiltonian, 



E n iip n i m (r) 



(5) 



H=2m+ — +V Q (r), 
m 



(6) 



is determined from combining two heavy quark propagators with each other. 

We wish to study relativistic corrections to the v = limit; after a Foldy-Wouthuysen- 
Tani transformation, the Feynman propagator is expanded in terms of the heavy quark 



Everything is consistently rotated to Euclidean space-time. 
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velocity^, v, around the static solution, in order to determine the propagator K(x, y) = 
q(x)q'(y). To 0(v 2 ) the propagator is given by, 

— d±K{x, y) = {igA^ + m — 5m(fi, m) + 0(x)) K(x, y), (7) 

with the well known terms, 

D 2 4 

= ~ V^Ci^.mjOi^m), (8) 

i=l 

1 fom) = ^j r , (9) 
2 (n;m) = ^-a-B, (10) 

3 (/x;m) = -i|M(D-E-E-D), (11) 
orrr 

4 (/i;m) = ^| ff .(DxE-ExD). (12) 

and 5j are color-electric and -magnetic field components. The heavy quark two-spinor 
q(x) consists of the large components of the original Dirac four-spinor after the Foldy- 
Wouthuysen-Tani rotation. 

Since we have truncated the expansion at a fixed power of v, we have lost renormalizabil- 
ity and the ultraviolet behavior is changed with respect to QCDf]. The theory is only effective 
and valid in the range of small gluon momenta q < \i. Whenever the Oj are determined at 
a scale \i that differs from m, the couplings m) (that are unity at tree-level) have to be 
adjusted by matching the effective theory to QCD at this scale; this guarantees the condition 
Ci(m,m) = 1 to hold. The zero point energy is shifted by 5m(fi,m) = E (fi) — E (m), with 
respect to QCD, where the static quark self-energy Eo(p) can be estimated from perturba- 
tion theory [PU|-[2^|. Due to this self-energy, the pole mass is shifted in respect torn- 5m 
within the propagator: m pole = m — 5m + E (fi) = m + E (m). 

Note that the Hamiltonian which corresponds to K is identical to that of NRQCD to 
order v 2 (up to irrelevant terms that are introduced to remove doublers and stabilize the 
evolution of the propagator on a discrete lattice with spacing a). In the case of NRQCD, the 
effective lattice theory is matched to continuum QCD in one step, such that the coefficients 



2 Formally, this procedure is equivalent to expanding the Dirac equation in powers of 1/c where c 
denotes the speed of light. Note, that in some of the NRQCD literature our 0(v 2 ) corrections are 
counted as 0(v 4 ). 

3 This fact gave rise to a discussion on a supposed discrepancy between the Eichten-Feinberg- 
Gromes results |7|,^] and perturbative expansions |23-p5|1 in powers of the coupling, g, where addi- 



tional terms that depend logarithmically on the mass occur after regulating loop diagrams. These 
terms are now understood to arise from changes in the ordering of integrations, and the underlying 



problem is resolved [1C]. 
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Ci do not only depend on m/n (or ma) but also on the lattice coupling g(a). We start 
from an effective theory, formulated in the continuum, such that the matching procedures 
(QCD-effective theory and continuum-lattice) will be treated in two separate steps. 

In addition, the operator, acting on K in Eq. (0), has the same structure to order 1/m 2 as 
the Lagrangian of heavy quark effective theory (HQET). Therefore, the matching coefficients 
can be taken from Refs. |25H3l||Fl 



c 2 (//,m) = , (13) 




2r, 



\a s {m) J 

C4(//, m) — 2c 2 (/i, m) — 1. (15) 

In order to evaluate masses of heavy quarkonia, we have to combine a propagator of 
a quark qi of mass mi with one of an antiquark q 2 of mass m 2 . In following the steps of 
Refs. |9|,|1(J, one can obtain the nonrelativistic Schrodinger-Pauli Hamiltonianf] (in the cm. 
system, i.e. p = Pi = — P2 and L = L 1 = L 2 ), 

H = V m; - 8mi + c^m;)— -3 + V{r, p, L, Si, S 2 ), (16) 

i= i \ Zm i % m i ) 

where the potential, 

V(r, p, L, Si, S 2 ) = V(r) + V sd (r, L, S 1; S 2 ) + V md (r, p), (17) 

consists of a central part, SD and MD corrections. 

Note that under renormalization group transformations the spin-spin interaction term 
of the effective two-particle Lagrangian [Jd^x (q\a ■ Y$qi)(q\<J ■ Bg 2 )] undergoes mixing with 
two local dimension six color singlet two-fermion terms that have to be included at order 
1/m 2 into the effective heavy quark Lagrangian. This very fact gives rise to the radiative 
correction term pllfl, 



In Refs. |26, 31,p2[, it has been shown that the kinetic energy term — D /(2m) does not undergo 
renormalization. Unlike in lattice NRQCD, where the quark mass becomes multiplicatively renor- 
malized [27|, here the mass does not enter as a dynamical variable of the simulation, but rather as 



an expansion parameter. The correction to the kinetic energy, that contains the dimension seven 
operator q^~D 4 q, however, is accompanied by a nontrivial coefficient c\(fi, m). This coefficient as 
well as the mixing matrix between Oi(fi;m) and lower dimensional operators has not yet been 
determined. For this reason, for the time being, we assume c\ ~ 1. 

5 The derivation of this expression from QCD is nontrivial ||. 

6 We have substituted the factor 87rC^a s (^t)5 3 (r) of the reference by the potential V^r), which is 
equivalent at this order in a s . 
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1 3 

SV sd = -^^SiS 2 -c 2 (m 2 , mi) (l - c 2 (/i, m 2 )) Vk(r), (18) 

within the SD potential below. 

The complete result on the potential to order v 2 with one-loop matching coefficients 
turns out to be, 

V(r) = V (r) + £ -^c 3 (m,) (vV (r) + VV a E (r 

i=l ° m i 

-E^NW^r), (19) 
i=i 

V«r, L, Sl , S 2 ) = + * ) L g!tlg™ 
, Si + S 2 l c + ^(r) 

171x1712 r 

+ ^L C2 (m 1 )c 2 (m 2 ) J R i ,\/ 3 (r) (20) 
mim 2 

+ ^2(m 1 )c 2 (m 2 ) - ^c 2 (m 2 ,mi) 1 - c 2 2 {m 2 ) ^ \/ 4 (r) 

+ / Si _ S2 \ L c-[V{(r) + V{(r)] 

\rrii 171% j r 
| Sx- Sa L c-^(r) 
mii7i2 r 



and 



Knd(r, p) = — { M , [ W) - i^K(r)]} (21) 

T(l\T(l 2 



2 1 



Weyl 



with 



= ^ " |. (22) 



1 

c± = c ± (ii,m 1 ,m 2 ) = - [c 2 (ju,rai) ± c 2 (/-t,m 2 )] , m x > m 2 , (23) 

Cj(m) = Ci(fi,m). (24) 

The symbol {a, 6, c}w cy i = ^{a, {&, c}} denotes Weyl ordering of the three arguments. 
V{, . . . , V4 are related to spin-orbit and spin-spin interactions. The MD potential gives 
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rise to correction terms of the form -L 2 , -4L 2 , -p 2 , - and 5 3 (r). The correction to the static 
potential includes, besides V 2 V^ and V 2 V^ , the expected Darwin term V 2 V . Note that 
V{, . . . , V4 as well as V 2 V^ and V 2 V^ depend on the matching scale /1 while Vo as well as 
Vb, ■ ■ ■ , V e are scale independent^]. In what follows, we will refer to the functions V[, . . . , V4 
as SD potentials, V&, . . . , V e as MD potentials, and V 2 V^ and V 2 V^ as corrections to the 
central potential. 

In order to derive the Hamiltonian from one-particle propagators, one has to assume that 
interactions between the two quarks are functions of a single global time coordinate (instan- 
taneous approximation). Unlike NRQCD, the above Hamiltonian cannot be generalized to 
higher orders in v since this would involve higher than first order temporal derivatives of 
the quark momenta, which, on the quantum level, cannot be reexpressed in terms of the 
canonical coordinates. 

Vo, V 2 Va , V 2 V^ ', V{, . . . , V4 and Vb, ■ ■ ■ , V e can be computed from lattice correlation 
functions (in Euclidean time) of Wilson loop like operators. Due to Lorentz invariance, 



certain pairs of potentials are related to the static potential by the Gromes [g3| and BBMP 
relations, 



V£(Ai; r) - V/Gu; r) 
V b (r) + 2V d (r) 



^o'(r) - \v {r), 



V c (r) + 2V e (r) = --V '(r) 



(25) 
(26) 

(27) 



such that three potentials, e.g. V{, Vd and V e can be eliminated from the Hamiltonian. 
From arguments, similar to those of Ref. ||, it is evident that the combination V^(fi;r) + 
2V 2 V^ B (/i;r) is a function of the static potential and thus scale independent. Given this 



observation, the structure of the Hamiltonian [Eqs. (|iq)-(]20|)l and the Gromes relation 
[Eq. (|25|)1, we can deduce the following one-loop relations between potentials, evaluated at 
cut-off scales \i\ and /x 2 , 



W/G^r) 



c 3 (/ii,/i 2 )V 2 Kf (/x x ;r) + ^(/ii,^) - 1] V 2 V (r) 
V 2 Vf( /Ul ;r) + ^V 4 ( / i 1 ;r) 



C 2 (/^l,^2) 



W tt B ( W ;r)- 
v?an;0-[i 

c 2 (/ii,/i2)U 2 '(^i;r 
03(^1, ^2)^3(^1; r 



W a fl (/* 2 ;r 
V^(/x 2 ; r 

^3(^2; r 

V4(// 2 ;r) = - 7c 2 (/ii,^ 2 ) - 3 V 4 (//i;r). 



1 - c 2 2 (fn, fi 2 )\ VtQwr), 
c 2 (^i,/i2)] V£(fii;r), 



(28) 

(29) 

(30) 
(31) 
(32) 

(33) 



7 The latter potentials originate from perturbing a quark world line, along which the field A4 of 
Eq. (|?]) contributes to the propagator, around the classical particle trajectory. Since an overall 
renormalization of the gluon fields can be absorbed into the quark wave function normalization, 
Vb, ■ ■ ■ ,V e are scale independent (like Vb). 
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B. How to compute the potentials 



In the Schrddinger-Pauli approach, introduced above, the quarks interact through a po- 
tential that only depends on the distance, spins, and momenta of the sources, Eqs. ( |17D - 
(PI). The time dependence has been separated and is implicitly included into coefficient 
functions of various interaction terms, the central, SD and MD potentials. These can be 
computed by a nonperturbative integration over gluonic interactions. Therefore, the po- 
tentials incorporate a summation over all possible interaction times, t. One obtains the 
following expressions in terms of expectation values in presence of a gauge field background 
for the corrections to the static potential H] 8 , 



c 



V 2 Vf (R) = 2 lim f di«E(O,0)E((M)» 

J 

V 2 Vf(R) = 2 lim [ T dt ((B(O,0)B(O,t)))vK 



where the superscript "c" denotes the connected part, 

((£i(ni, 0)4-(n 2 ,t)»^ = (^(in.O)^!!^)})^ 

- lim((£ i (n 1 ,0)4-(n 2) t , )>>w- 

t — >oo 



For the SD potentials one finds fTJH) 

2e ijk lim fdt t ( (Bi (0, 0) % (0, t ) ) ) w , 
e iifc lim / T dt t ( (A (0, 0) Ej (R, t) ) ) w , 



f Vi'(R) 



^v 2 '(r) 



RijV 3 {R) 



2 lim 

r— »oo 



dt 



((A(0,0)4(R,t))) w 



«B(0,0)B(R,t)» 



if 



V 4 (R) = 2 Jim / dt((B(0,0)B(R,t))) w . 
Finally, the MD potentials are 

Vfc(R) = - 



1 



lim / dtt 2 ((E(0,0)E(R,t))) c w , 

RijVciR) = lim fdtt 2 [((^(0,0)4(R,t))) c w 
Si 



^«E(0,0)E(R,t)» 



(34) 
(35) 



(36) 

(37) 
(38) 
(39) 

(40) 

(41) 
(42) 



We have recast all expressions into forms that are more suitable for lattice simulations. Via 
spectral decompositions of the underlying correlation functions, equality between our definitions 
and those of Refs. (7|-|9| can easily be shown. 
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V d (R) = - hm / dtt 2 ((E(0,0)E(0 } t))) c w , 



— lim 

5iA 



dtt 1 



((E i (O,0)E j (O,t))) 



c 
U 



(43) 
(44) 



((E(O,0)E(O,t)» c w 



R denotes a lattice vector of length R = ra 1 . At small lattice spacing a, the above potentials 
should approach their continuum counterparts and rotational invariance is expected to be 
restored, V (R) = aV (r), V( >2 (R) = a 2 V{ 2 {^r), V 3A (R) = a*V 3A (fi;r), V 2 V a E ' B (R) = 
a 3 V 2 V a E > B (fi] r), H, c , rf , e (R) = aV b>cAe (r), where // = n/a. 

Throughout the previous equations, the expectation value ((FiF 2 )) w is defined as, 

(TrP[exp {igJ dw dx tl A tl )F 1 F 2 }} 



({FiF 2 )) 



w 



(Tr V [exp (ig J 9W dx^ A^)}) 



(45) 



where dW represents a closed path [the contour of a Wilson loop VT(R, T)\ and V denotes 
path ordering of the arguments. Although we have chosen a lattice inspired notation for 
the potentials [Eqs. (p4[)-(^)], so far everything is generally applicable to lattice as well as 
continuum formulations of QCD. In following Huntley and Michael (HM) |)TB| , we implement 
the discretized version of Eq. (f45|), 



((FiF 2 )) 



w 



(V 


w (nx-n^na-nl)^ 


)(W) 


(v[ 


+ nj)" 


) (v [w(u 2 + nj)" 


> 



(46) 



where the subscript i = 1,2 represents the multi-index (n^/ij,^) and are integer valued 
four-vectors. The subscript "tl" indicates that only the traceless part is to be taken: (A) t i = 
A — |TrA Fi are related to the electric and magnetic fields in the following way, 



\J,V 



ga 2 Fn U , Ei 



i4, 



T&ijkFjk- 



(47) 



These conventions eliminate imaginary phases and factors g 2 a from Eqs. (|34])-( 

We have taken Ii^ v {n) to be the spatial average of the four (two) plaquettes, enclosing 
the lattice point n for magnetic (electric) fields, 



and 



with 



ny(n) = - [Pi,j(n) + Pi,-A n ) + P-i,-j( n ) + p -i,A n )\ 



n i4 in + -4 = - [P iA (n] + P_ M (n)] 



(48) 



(49) 



P^{n) = U^n)U u (n + fi)UUn + v)Ul(n), U.„{n) = UUn - /}). 



(50) 
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With this choice of II, Eq. ([46]) is correct up to 0(a 2 ), the discretization error of the Wilson 
action, used for generating the gauge field background. Note that the electric fields are 
living at half-integer time coordinates, in between two adjacent spatial lattice hyperplanes. 
U^n) is the SU(3) link variable, related to the field A^x) at x = (n + ~/t)a: U t 



n 



Pexp igaf™ +t *dn' u A u (n'a) 

In practical computation, the temporal extent T of the Wilson loop W [within Eqs. 
([44])] is adapted according to the formula T = t + At 1 + At 2 - Aij, the separations of the 
"ears" F\ and F 2 from the corresponding spatial closures of the Wilson loop, are kept fixed 
throughout the simulation while the interaction time t, is varied. The discretized version of 



the nominator of the correlation function Eq. (45) is visualized for the case of an electric 
and a magnetic ear in Fig. [I]. Strictly speaking, Eqs. (J53)-(K1) apply in the limits Ati — > oo 
only. Ati (Ai 2 ) represents the time the gluon field has to decay into the ground-state, after 
(before) creation (annihilation) of the qq^ state and is a control parameter of the simulation. 



C. Theoretical expectations 

1. General considerations 

In addition to the exact Gromes and BBMP constraints, Eqs. (p5])-(p7|), some approx- 
imate relations between the SD potentials are anticipated from exchange symmetry argu- 
ments. We start from the standard assumption that the origin of the static potential is due 
to vector- and scalar-like gluon exchange contributions. Given that a vector-like exchange 



can grow at most logarithmically with r fl34| , the nature of the linear part of the confining 
potential can only be scalar. As we will see, V^'(r) is short ranged, such that the confining 
part only contributes to V{(r). This leads us to expect V%(r) to be purely vector-like. Under 
the additional assumptions that pseudoscalar contributions can be neglected and that V[ 



does not contain a vector-like piece, one ends up with the scenario of interrelations |33 



V 3 (r) = ^ - V?(r), (51) 
V 4 (r) =2VV 2 (r), (52) 

which of course has to be in agreement with leading order perturbation theory. However, 
Eqs. fl5l|) -(|52"D hold true for any effective gluon propagator that transforms like a Lorentz 
vector. Unlike the Gromes relation, the above relations cannot be exact, which is evident 
from Eqs. (BT|)-(B3). 



2. One gluon exchange potentials 

In order to parameterize the short range behavior of the potentials, it is useful to resort to 
weak coupling perturbation theory. For modeling of lattice artefacts, we have calculated the 
SD potentials to tree- level in Ref. [^TJ. Here, we supplement these results by the remaining 
0(v 2 ) potentials. 

Lattice potentials. In the following we will use the conventions, 
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G L (R) 



LI 



E 

o q^O 



3 iqR 



2 ' 



cji = 2 sin 



and 



with 



*i(R) = f!fE ' 



2tt 



-mi, rrii 



c q^O (Si Qi ) 



2\2> 



+ 1 — 



(53) 



(54) 



(55) 



L a denotes the number of lattice sites along a linear spatial extent. Note that the above 
functions have the large- R behavior (for L a ^> R), 



G L (R) - ^, F L (R) - F L (0) - 

where -Fl(O) diverges linearly with 
We find 

V (R) = -C F g 2 (G L (R) - G L (0)) , 



V 3 (R) 



R 2 



RjR, 



-C F g 2 A J A k E l G L (R) 



(56) 

(57) 
(58) 

(59) 



j -fife 



for jj^k,ij^j,ij^k and i?j ^ 0, R k ^ 0. Unless explicitly stated, no summations over 
indices that appear twice are performed. For 2i? 2 ^ R 2 + R\, k as above, we can derive 
the expression, 



V, = 



R 2 



o p2 d2 p2 



^I[4Af)S-Af(S fe + S)-Af(S, + S) 



The remaining potentials are given by, 

^ 4 (R) = -2C F /£A?^ G L (R), 

C F g 2 



V b (R) = ^A{ 2) S i F L (R) + 6~G L (R)) , 

K(R) = 



(60) 

(61) 
(62) 



i? 2 



2 ZIURj - 5ijR 2 



-3A,A 7 F L (R) 



+ Sij A {2) ~ k F L (R) + 6 (S - SO G L (R)) 



(63) 



Vd(R) = 



4 



Gl(0) + Gl(1) + \f l {2) - l -F L {Q) 
G L (0) {L a -> oo), 



V 2 Vf (R) = -3C F # 2 [G L (0) - G L (2)] « -0.629525C F £ 2 , 
V 2 Vf (R) = -6C F g 2 [G L (0) - G L (V2)} « -1.185237C F ^ 2 



(64) 

(65) 
(66) 
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with f(l) = f(l), f(y/2) = f(l + 2), f{2) = /(21). The numerical values refer to the 
infinite volume limit. In this limit, one obtains Gl(0) ~ 0.2527310. The potentials V{ and 
V e vanish to lowest order perturbation theory. 

The Casimir factor of SU(3) gauge theory is Cf = 4/3. The As and Ss denote the finite 
difference and averaging operators, 



Ai/(n) 
AfV(n) 

A (2) 

Si/(n) 



I[/(n + S)-/(n-l 
/(n + i)-2/(n) + /(n-i) ) 

i 

I[/(n + i)+2/(n) + /(n-i) 



(67) 

(68) 
(69) 

(70) 
(71) 



In Ref. |[2rj| , we have proven that an exact lattice analogue to the Gromes relation does 
not exist. However, the Gromes as well as the BBMP relations will be retrieved in the 
continuum limit and approximately hold within the scaling region on the lattice for R ^> 1. 

Continuum potentials. In continuum perturbation theory, one obtains the following 
tree-level expressions, 



V (r) 
Vj(r) 
V 3 {r) 
VJr) 



-C F a s 



dq 



3 piqr 



2tt 2 q 2 



dq 3 



—Cpd 

CfOI 



27r 2 q 2 r 
dq 3 (qr) : 
2tt 2 q 2 r 2 
dq 3 



M_ e iqr 



7T JO 

2a. 



3 iqr 



-c 



-a 



IT JO 

2a 



sin qr 
qr 

DO 

dq q 2 rji(qr) 
dq q 2 j 2 (qr) 



D «qr 



2tt 2 



C, 



2a, 



7T 



dqq 



IT JO 

2 sin qr 



qr 



(72) 
(73) 
(74) 
(75) 



g 2 /(4iT). Vb and V c are given by — 2Vo/3 and Vq/2, respectively. The self- 
Cf5 ,2 Gl(0) has been subtracted from Vq. V[ and V e vanish to lowest 



with a s 

energy Cy = a^ 1 

order perturbation theory while V 2 Vq E , ^Vj 3 and only contain diverging self-energy 
contributions. A linear confining contribution can be introduced by adding a — l/g 4 -term to 
Vq in momentum space. The integrals for the SD potentials are suppressed in the infrared 
region like q 2 or g 3 , such that we naively expect perturbation theory to be more reliable in 
this case than for the static potential or Vb and V c . 
Eqs. ©-fD yield, 



V (r) = -Q 



a. 



a* 



V 3 (r) = 3C F ^, 



(76) 
(77) 
(78) 
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V b (r) 
VJr) 



8nC F a s 5 3 (r) 

3 r 
2 r 



(79) 
(80) 

(81) 



in agreement with the large- R (i.e. r>o) expectations of Eqs. 



[cf. Eq. ©]. 



3. Large distance behavior 

In combining the large-r behavior from the minimal area law (MAL) of fluctuating world 
sheets (including a perimeter term) with the expectation from tree-level perturbation 
theory, Eqs. QTEp-flglp, one obtains, 



V (r 
W a E (r 
V{(r 
V 3 (r 
V b (r 
VAr 



V c h «r, 

r 



V!(r) = ^. 
3^, \/ 4 (r) = 87r(e-/i)5 3 (r), 

26 /t -t r / \ \ €■ K 

Ch H r, = r, 

b 3r 9 ' V ; 2r 6 ' 

Cd-^r, V e (r) = --r, 



(82) 
(83) 
(84) 
(85) 
(86) 
(87) 



with e = Cpa a , in agreement with the Gromes and BBMP relations, Eqs. (P5])-(|27|). From 
the perimeter term in MAL one obtains = = C b = and Cd = —Cyj^. In tree-level 
perturbation theory, one finds a consistent infinite volume result Cd = —Cy/4: [Eq. (|64T)1, 
while comes out to be significantly smaller than C^ [Eqs. (j55|)-(]6"6D]. However, our 
numerical data shows ~ , in agreement with the MAL result. From Eqs. fl2"8"|)-(|3"0"D, 
it is obvious that the tree-level perturbative expectations Eqs. (|76|) -(^TD cannot adequately 
describe the potentials at all scales, fi. In general, V{ and V 2 ' will undergo mixing, such that 
V{ will attain a Coulomb-like contribution. For the same reason, V 2 V^ is expected to include 
a 1/r-piece. We have accounted for this fact by allowing for two additional parameters, b 
and h. In principle, V 2 V r Q E ' B can also contain 5-like admixtures [Eqs. (p8"D, (|29|)1, which we 
have ignored in Eq. fl83f) . 



III. LATTICE SIMULATIONS 



In Ref. |Tl|, we have developed suitable techniques for a lattice evaluation of the poten- 



tials and applied them to SU(2) gauge theory. We investigated possible sources of systematic 
errors such as finite size effects. In this section, we describe details of our SU(3) simulations, 
insofar these differ from the SU(2) study. 



14 



A. Simulation parameters 



We analyse two sets of Monte Carlo configurations that have been generated with the 
standard Wilson action on hypercubic lattices of volumes V = L\L T = 16 4 at (3 — 6.0 and 
V = 32 4 at f3 = 6.2 (Table |). The above couplings correspond to inverse lattice spacings 
a^ 1 pa 2.1 GeV and a -1 pa 2.9 GeV, respectively. The scale has been determined from the 
value ^[k = 468 MeV for the string tension that we obtain from the fit to the bottomonium 
spectrum of Sec. |VT]. The number of independent Monte Carlo configurations n con f, gener- 
ated at each set of parameters, is included into the table. Based on previous experience, we 
expect finite size effects to be below statistical accuracy at these volumes [ITTJJT^,pTJ . For the 



updating of the gauge fields, a hybrid of Fabricius-Haan heatbath |35| and an overrelaxation 
algorithm has been implemented ||36|| . Within both procedures, we successively update the 
three diagonal SU(2) subgroups of a given link p7| . The heatbath sweeps have been ran- 
domly mixed with overrelaxation steps with probability 1/7. The links have been visited in 
lexicographical ordering within hypercubes of 2 4 lattice sites, i.e., within each such hyper- 
cube, first all links pointing into direction 1 are visited site by site, then all links in direction 
2 etc.. After 2000 initial heatbath thermalization sweeps in either case, measurements are 
taken every 100 sweeps to ensure de-correlation. We find no evidence for any autocorrelation 
effects between these configurations. 

B. Noise reduction 

Statistical fluctuations have been reduced by "integrating out" temporal links that ap- 
pear within the Wilson loops and the electric ears analytically, wherever possible. By "link 
integration" we mean the following substitution [|38f . 

1 dZ L um dUUe s ^W 

U ' M - lUn) = zem = W» • (88) 

with 

S nifl (U) = Tr (F»C/t + UFl(n)) , (89) 

and 

F^n) = U u (n)U^(n + u)Ul(n + ft), (3 = \. (90) 

Wi{n) is in general no longer an SU(3) element. In this way, time-like links are replaced 
by the mean value they take in the neighborhood of the enclosing staples F^in). Only 
those links that do not share a common plaquette can be integrated independently, without 
changing expectation values. 

We have attempted to compute W±{n) analytically for the case of SU(3) gauge theory. 



Based on the character expansion of SU(N) matrices of Ref. |39|], the following expression 
can be obtainedPI, 



9 For simplicity, we suppress spatial coordinates and Dirac indices of U, F and W. 
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with 



Z(F) = J su3 dU exp (Tr (FU ] + [7F f )) 

dx f dy ( Fix) 

— f — exp [xQ + y + 

2m J 2m \ xy 



Q = det(F) + det(F t ), 
P(x) = 1 + Tr(FF f )x + - [Tr 2 (FF t ) - Tr{FF^) 2 } x 2 + det(FF f ) 



(91) 
(92) 

(93) 
(94) 



From 



/A nr> pQ% 



K n [0] = f±O^M2RM), 



(95) 
(96) 



where I n denote the modified Bessel functions and R{x) = JP(x)/x, one obtains |4C} 



Jo ' W - Jl m +K ° 



dP 

apt 



(97) 



such that 



W 



Jn 



J X G + K X F + K 2 F (Tr(F f F) - F ] f) + K 3 det(F)G 



(98) 



where K n = K n [l), Gu = ^€ijkei mn F* m F^ n . Note that J n and K n are real numbers. For the 
computation of Bessel functions we use the asymptotic expansion, 



In(z) 



V2 



2vrz p K ' z3 



(99) 



with 



Aj(n) 



{An 2 - l 2 )(4n 2 - 3 2 ) • • ■ {An 2 - (2j - l) 2 ) 



(100) 



up to fifth order in j. The above expansion is valid for arguments, z = 2R{x), with large 



modulus. A circular integration path with radius \x\ = 0.015 turns out to be appropriate |4T 



at (3 ~ 6. A Gaussian quadrature algorithm with 64 abscissas is used. By exploiting the 
symmetry of the contour integrals Eqs. (|95|) and (|9~6D under the transformation x — > —x, we 
are able to reduce the computational effort by a factor 2. 
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C. Ground state enhancement 



In this section, we will discuss the control of excited state contributions at finite de- 
excitation times Aij. We found At = 2 to be appropriate for magnetic ears and At = 3/2 for 
electric ears (see Fig. |l|). The spatial transporters within the Wilson loops have been smeared 
to suppress excited state pollutions from the very beginning. Our smearing procedure p2|j42 



consists of iteratively replacing each spatial link Ui(n) within the Wilson loop by a "fat" 
link, 

Ui{n) ]oiUi{ri) +J2 u j( n ) u i( n + j) U ji n + i) \ , (101) 

with free parameter a. M denotes an operator that projects the argument back into the 
SU(3) group: U = Af(A) E SU(3) with Re Tr{A^U} = max. Within this procedure, the 
(spatial) links are visited in the same lexicographical ordering as within the Monte Carlo 
updating of gauge configurations. We find satisfactory ground-state enhancement with the 
parameter choice r^ter = 100 and a = 2. 

From expectation values of Wilson loops, the static interquark potential can be deter- 
mined in the limit of large T, 

(W(R, T)> = C e^»T f 1 + j2 , (102) 

V n>0 °0 / 

where AV n = V n — V denotes the gap between the ground-state and the nth excited state 
(hybrid) potential. The R dependency has been omitted from the overlap coefficients < 
C n < 1 and potentials V n , AV n . The smearing procedure results in an increased weight C 
(with respect to C n , n > 1). In Fig. |2| the resulting static interquark potentials at (3 = 6.0 
and (3 = 6.2 are displayed. All ground-state overlaps turn out to be well above 0.8. 



Previous authors JT^|T^|T^j2y] have replaced the integrals over interaction times by dis- 
crete sums. This results in cut-off errors due to the finiteness of the integration bound r 
[Eqs. (|34|)-(P4D1 as well as 0(a 2 ) integration errors. Both sources of systematic uncertainties 
can be studied and reduced by exploiting transfer matrix techniques. In the following, we 
will briefly summarize some of the results we have already presented in Ref . [ETJ . The ratio 



of correlation functions between eared Wilson loops, Eq. (|^), is given by0, 

((FiF 2 ))w = E D^e-^ [i + E%e- A ^ M + •••]. (103) 

m 

All constants are understood to depend on R. For details on and E^, that are functions 



of the spatial positions and the color-electric or -magnetic components F\ and F 2 , see [^T 
The unwanted excited state contributions are suppressed by factors \E^\ < JC\/Cq as well 



10 The formula applies to the case At = At\ = At2- It remains valid on a qualitative level for 
Aii / At 2 with At = min{Ati, At 2 }. 
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as by 



-AViAi 



The smallest value of At that appears within an integral over interaction 



times will determine the reliability of the result. The bosonic string picture yields the large- 
R expectation AVi(R) = tt/R |43| for the lowest lying hybrid potential, which has been 
qualitatively confirmed in numerical studies |44],|45[ . 



The cylindrically symmetric creation operator that we use only projects onto states 
within the A\ g representation of the appropriate symmetry group @|- The lowest 
continuum angular momentum to which it couples is L = 0. The hybrid (L = 1) state E u is 
the next excitation [ II] . The operators used as magnetic ears have no overlap with the A lg 
state, such that all correlation functions that involve a magnetic ear decay exponentially 
with Euclidean time. This does not hold true for some of the correlators within the MD 
potentials and V 2 V^; those electric ears which are not orthogonal to R have a nonvanishing 
overlap with Ai g , such that the disconnected part in Eq. ( |103| ), D^ 2 , does not vanish and 
has to be explicitly subtracted in Eqs. and (|4l| ) - (|44]) . 

Note that the are not normalized and can be negative. However, due to invariance 
under time inversion, the correlation functions for V[ and [Eqs. ([37] and (0)] have to 
vanish at t — 0, such that J2 m = in this case. In combining Eq. ( |103| ) with Eqs. 
(H) or (H)-(gO]), we obtain, 



V 2 K 



E,B 



V 3 a oc 



m>0 



dt D 



12, 



-AV m t 



E 

m>0 



D 



12 



AV„. 



(104) 



with appropriate color field positions ni,n 2 and components /ii, v\, {12, ^2- Eqs. 
©yield, 

% oc E rdttD^e-^ = £ -S-j, 

m>0 J0 m>0 (AV^J 

while from Eqs. (p|)-(|4^) we obtain, 



V aM , d <x / dtt 2 D%e 



00 . ni2 

2 r,12-AV m t ^ ^ 



and 



(105) 



(106) 



The parameters Dj^ and AV m can be fixed from fits of the data to Eq. (|103|) . The 
hybrid potentials V m can in principle also be determined independently f44| . We leave this 
for future high precision studies on anisotropic lattices. For the time being, we evaluate 
the integrals Eqs. ([34])-(|44]) numerically. The interpolation method used for this purpose is 
inspired by the multi-exponential result of the spectral decomposition, Eq. 



D. Integration errors 

The 0(v 2 ) potentials are extracted from integrals over correlation functions [see 
Eqs. (|3^)-(f44"Dl that depend on the interaction time t in a multi-exponential way. In the 
following, Ciit) will denote the two-point function which has to be integrated out in order 
to determine a potential Vi at a given value of R. For i = 1,2, Cj(t) will be weighted by 
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an additional factor t [Eqs. ©-©], for i = b,c,d,e by t 2 [Eqs. ©-flU)]. Two different 
methods of interpolating Ci(t) in between the discrete t- values have been adopted: 

1. We perform local exponential interpolations, which are expected to yield the most 
reliable results: 



C t {t) 
Q(t + 1) 



(107) 



for t < t' < t + 1 and Ci(t)Ci(t + 1) > 0. Due to the multi-exponential character of 
the correlation function (or statistical fluctuations) the sign might change within the given 
interval. Thus, for Ci(t)Ci(t + 1) < 0, we interpolate linearly, 



Q(t') = d{t) + [d(t + 1) - Q(t)}(t' - t). 



(108) 



For Ci(t) and C 2 (t) quadratical interpolations are performed within the interval < t' < ~ 
to account for Ci(0) = C 2 (0) = 0, where we demand continuity of the interpolating function 
and its derivative at t — ~. 

2. In order to estimate systematic integration errors, simple linear interpolations of the 
data have been performed. 

All statistical errors have been bootstrapped. For each potential Vi, numerical integration 
has been performed up to a value t = Tj, with T{ chosen such that the result is stable 
(within statistical accuracy) under the replacement Tj — > Tj — 1 for all R. Systematic cut- 
off errors have been estimated from the exponential tail of fits to large-t data points and 
have been found to be negligible in all cases when compared to the statistical error from 
the numerical integration. Typically, Ti came out to be 6-8 lattice units For some of 
the correlation functions, the disconnected part had to be subtracted. Its value has been 
estimated by averaging data points within a range of t- values, t > rf, under variation of rf 
until a plateau was reached. Subsequently, the resulting value has been subtracted from the 
correlator, before proceeding to interpolation methods 1 and 2. In all cases, rf — 1 < Tj < r" 
was found, in support of self-consistency of the method. 

In what follows, we always state the result from the exponential interpolation method 
with the bootstrapped statistical error and a systematic error that corresponds to the differ- 
ence between the results obtained from the two methods. We find the systematic error to be 
the dominant source of uncertainty, which can only be reduced by decreasing the temporal 
lattice spacing. 



IV. THE MATCHING PROBLEM 
A. Renormalization of lattice operators 

The relativistic corrections to the static potential are computed from amplitudes of 
correlation functions rather than from eigenvalues of the transfer matrix. Therefore, they 
undergo renormalization. This is in accord with the fact that the electric and magnetic ears 
explicitly depend on the lattice scale a and, thus, discretization. 

As in the low energy regime of interest the renormalization constants are likely to receive 
relevant corrections beyond one-loop perturbation theory, we apply the nonperturbative HM 
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renormalization prescription |T8] [cf. Eq. The HM procedure is similar to the mean 

field inspired tadpole improvement program, advocated in Ref. [47]. However, instead of 



just dividing correlators of eared Wilson loops by the square of the average plaquette, 



n) 



(109) 



a more sophisticated combination is chosen; the various orientations of ears are taken into 
account, such that the remaining renormalization constants will only differ from identity on 
a three-loop [1 + 0(g 6 )] or two-loop [1 + 0(g 4 )} level for operators involving magnetic and 
electric ears, respectively. Details are discussed in Ref. . 



In the case of tadpole improvement, each electric or magnetic field gE or gB, appearing 
within the correlators of Eqs. (34)— (44), is multiplied by a constant tadpole = l/Uu. The 
difference between this procedure and the HM scheme can be parameterized in terms of a 
ratio, 



Q 



V 



2 (W) Uc 



(110) 



that depends on the orientation of the well as R and t. In Fig. [3] we display this 

ratio for all independent color-electric and -magnetic components at R = 4 and (3 = 6.2 
as a function of t. No significant t-dependence is observed, such that the renormalization 
constants within integrals Eqs. (33)-(f|lJ) factorizeQ all Q-factors saturate into asymptotic 
values for t > 4, within statistical errors of (9(10~ 4 ). In Fig. |], the R dependence of Q 
is depicted for on-axis separations. In the case of magnetic ears, the result appears to be 
rather insensitive to the component and R. For electric ears, Q changes significantly with R 
as well as the component. However, at large separations, the electric components approach 
a common value too. We call these plateau values Qb and Qe, respectively. In Table ||, 
we compare Z tadpo i e with the HM constants Z B = Z tadpo i e /Q B and Z E = Z taApoXe /Q E for the 
two (3- values. The magnetic renormalization constant differs only by less than 1 % from the 
tadpole value while for the electric fields this difference amounts to 3-4 %. As expected, 
the disagreement decreases with the lattice spacing a. We find it interesting to notice that 



the factors Q are smaller than ratios of nonperturbative values [53] of the coefficient of the 



clover term within the Sheikholeslami-Wohlert fermion action |49| and its tadpole guesses. 
Also, the correction goes opposite in the present case; Z B and Z B come out to be smaller 
than the tadpole estimate ^tadpole- 
Direct numerical checks of the accuracy of the HM approach are possible in two ways, 
namely (i) by varying the lattice resolution a and a scaling test of the results^ and (ii) by 



11 This behavior is expected from the spectral decomposition of the correlation function [21], 
where the renormalization factor corresponds to a constant, l/^oo' that only depends on R and 
the specifications of the ear i. Any residual time dependence has to be attributed to the fmiteness 
of At,-. 
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Due to the running of the matching constants to the full theory with the lattice scale, residual 
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comparing the data with predictions from the exact Gromes and BBMP relations, Eqs. fl25|)- 
(P7|)) between SD or MD potentials and the static potential (which does not undergo renor- 
malization) , 

^'renWa; r) ~ Vi^Tr/a; r) = V^t), (111) 

H, ren (r) + 2V d , Iea (r) = T -V^r) - \v*{r), (112) 

6 2 

V C)TW {r) + 2V eilea {r) = ~V^{r). (113) 

In Fig. [| we check our data on V{ — V{ against the force, obtained from fits to the static 
potential, Vo(r) according to the parametrization Eq. ( |117|) below. As can be seen, the two 



data sets scale onto each other and reproduce the static force. The BBMP relations are only 
satisfied on a qualitative level as Figs. |] and [7| demonstrate; substantial lattice artefacts are 
responsible for deviations from the expectations in the region of small R. 



B. Matching constants 

In order to calculate the matching constants between the effective nonrelativistic Hamil- 
tonian of Eqs. (|i~6|)-(pT|) and QCD, we require values for the strong coupling constant at 
scales q = it /a and q = mb,m c in a given renormalization scheme. We decide to use the 
"V" scheme of Ref . [pO] , and compute the running coupling from the average plaquette as 



suggested in Ref. [47 



Cl 



InC/r 







(-) 


+ 0.1058 


\aqj 





;ii4) 



where C\ = 1/3 and b = ll/(167r 2 ) for SU(3) gauge theory. 

We use the plaquette values of Table [TT] and Eqs. fllBp-flTop to obtain the matching 



constants c 2 (/i, m) and c 3 (/i, m) at (3 — 6.0 and (3 = 6.2, listed in Table [nT[ We have assumed 
rrib = 4.7 GeV and m c = 1.3 GeV for the bottom and charm quark masses, respectively. The 
value y/H = 468 MeV has been used to fix the lattice scale. These values are obtained from 
the quarkonia spectroscopy below. Since ay(q) depends only logarithmically on q, accurate 
values for a, m c and m b are not required. 

In the case of bottomonium all constants turn out to be reasonably close to 1, such that 
one-loop perturbation theory appears to be trustworthy. However, for m = m c , the size of 
the constants indicates that higher order corrections cannot be neglected at present lattice 
spacings a 7r/m c . Increasing the lattice spacing would result in increased lattice artefacts 
as well as a larger uncertainty in the renormalization factors Zb and Ze that relate the 
lattice potentials to their continuum counterparts. In order to achieve a reasonable balance 
between the uncertainties involved in both matching procedures for the charmonium family, 
improved lattice actions [BTJ would have to be considered. 



scaling violations for the SD potentials are expected from Eqs. ( |30[ ) and (|3l|). 
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V. RESULTS ON THE POTENTIALS 



We present numerical results and parametrizations on the static potential, the relativistic 
corrections to the central potential, and the SD and MD potentials. We compare the short 
range SD potentials V 2 ', V 3 and V4 and the MD potential V c to lattice perturbation theory. 
The short range SD potentials might provide another access to the determination of the 
QCD coupling ay(ir/a), quite in the spirit of the role of the fine structure constant in the 
analysis of atomic level splittings. 



A. The static potential 

The lattice potential Vo(R) has been computed from smeared Wilson loops by use of the 
method described in Ref. f22]| . Our general strategy is to derive interpolating parametriza- 
tions of the lattice data points which will enable us to compare the results to continuum 
expectations. Weak coupling continuum and lattice predictions on the potentials have been 
presented in Sec. [II] [Eqs. (|76|) -(|8lD and Eqs. (|57|)-(|66D, respectively], such that we can cor- 
rect the lattice data for the differences between both tree-level forms before attempting to 
fit them to a continuous parametrization. Let 

V ,cont(R) = VoQt) - gSV (R) - Cy, (115) 

with 

<5Vo(R) = -47rG L (R) + 4 (116) 



be the tree- level corrected static potential. Gl(R) is the lattice gluon propagator of Eq. 
The static lattice potential is fitted to the five-parameter ansatz [including g and Cy of 
Eq. ( |115| ) as fit parameters], 

V 0>cwt (R) = KR-^ + ^ (117) 

with string tension K = na 2 and Coulomb coefficient e. The 1/R 2 correction, that accounts 
for the running of the coupling, is not meant to be physical but has been introduced to 
effectively parameterize the data within the given range of r-values. The resulting parameter 
values are displayed in Table IV ^. For technical reasons related to the link integration 
procedure, only potential values for R > \/2 have been obtained, such that the fits do not 
include R = 1. 

In Fig. H, the potential Vo iCO nt from both /3-values is displayed in physical units (as 
obtained from y/~K = 468 MeV), together with a fit curve that corresponds to the (averaged) 
values of fit parameters e = 0.321(6) and f = af = 0.0082 (8) /\/k. As can be seen, the two 



13 The reduced x 2 ~ va l ues stated in the table do not take account of correlations between data 
points obtained at different R. 
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data sets scale nicely onto each other. Violations of rotational invariance are removed by 
the correction method, even at very small values of R, and the data is well described by the 
parametrization over the entire r-range. 



B. Corrections to the central potential 

Fits of V 2 V^ to the parametrization of Eq. 



V a V? = 3?--^ (118) 

with parameters and b have been performed. The resulting potential V 2 V r Q E — Cf is shown 
in Fig. |8|. Results on b are displayed in the first row of Table [V[ Systematic errors from the 
integration procedure are included in square brackets. We find the values \ 2 I^df = 0.5 
and x 2 /Ndf = 2.0 at (3 — 6.0 and (3 = 6.2, respectively, for the fit range R > y/2. These x 2 - 
values refer only to the statistical errors. The fitted curve that corresponds to the averaged 
value b = (0.86 ± 0.05 GeV) 2 is included into the figure. From Eq. (p8[) and the matching 



constants of Table pj, we expect scaling violations of about 10 % between the two data sets. 
Apart from the region r < 0.15 fm, which is polluted by lattice artefacts, this effect cannot 
be resolved within statistical accuracy. 

In Fig. [| we display V 2 V^ in lattice units at (3 = 6.2, where the error bars of this 
plot refer to the statistical uncertainty only. The (3 = 6.0 data exhibits the same quali- 
tative behavior. The large-i? data can be parameterized by a constant. Deviations from 
this constant at small i?-values, which are hidden within the systematic uncertainty of the 
integration, can be due either to lattice artefacts or to a tiny 5-like admixture that one 
might expect from Eq. (p9|). The numerical values (with statistical and systematic errors) 
are (7f = -1.02(1) (27) and <7f = -0.93(1) (24) at (3 = 6.0 and (3 = 6.2, respectively. These 
values have to be related to Cf = -1.00(2)(8) and 6? = -0.92(1)(8), such that Cf = Cf 
within errors. 

We conclude that the corrections to the central potential agree reasonably well with 
the expectations of Eq. ( [S3]) , with a parameter b w 4k w 0.9 2 GeV 2 . The strength of the 
effective Coulomb coupling is increased by about 2 % in case of the T family and 35-40 % 
for J /if) states, due to these correction terms. The self-energy type contributions to V 2 V^ 
and V 2 Vj 3 cancel each other at the present level of statistical accuracy. 



C. Spin-dependent potentials 



Our results on the first spin-orbit potential V[ are displayed in Fig. [10]. The two data 
sets show approximate scaling behavior. In addition to a constant long range contribution, 
—K, we find an attractive short range contribution that can be fitted to the Coulomb-like 
ansatz of Eq. (R3|) , 

V{(r) = -±-K, (119) 
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in agreement with our SU(2) investigation For these one-parameter fits we have con- 
strained the constant long range part to the value of the string tension, as obtained from 
the static potential. We find the values h = 0.071(12) and h = 0.065(9) for (3 = 6.0 and 
(3 = 6.2, respectively. As expected from Eq. fl30"D , h tends to decrease with (5. 

Taking the Gromes relation and the running coupling improved effective parametrization 
of Eq. (|117|) into account, we expect, 



Note, that we have added a l/r 3 -term to the expectation Eq. (|84| ) which accounts for a 
weakening of the effective coupling with decreasing source separation. From Eq. flBTD we 
expect the parametrization, 

m = M _ v (121) 

to approximate V 3 . 

Prior to comparing the data to the above continuum parametrizations, we attempt to 
correct for lattice artefacts. For this purpose we define, 

47T - 1 

^ 2 '(R) = ^-^^(R) " W (122) 
5V 3 (R) = J^t/ 3 ,t r ee(R) - j^. (123) 



The lattice tree-level potentials V^ tree and V3 :tree are defined in Eqs. (|5q)-(|60D. We then 
correct for lattice artefacts, 

^'cont(^) = " fla^CR), (124) 

V 3 , cont (R) = V 3 (R) - g 3 5V 3 (R), (125) 

and fit the potentials to the following ansatze: 

vl^(R) = ( 126 ) 

V 3 , cont (R) = 3 -§-J^> (127) 



where <?j, q and /, are fit parameters. The resulting parameter values are shown in Table [VI . 
Again, the x 2 - v alues refer to the statistical errors only. 

The fitted values f% and ^3 are in agreement with / as extracted from the static potential. 
Also, C2 and c 3 agree with e — h as computed from Vq and V\ reasonably well. Only the 
coefficients of the correction terms, gi and g 3 , come out to be about a factor 2 smaller 
than in the case of the static potential. The spin-orbit potential V^ cont and the spin-spin 
potential V 3tCont are displayed in Figs. [11] and [12], respectively, together with the theoretical 
expectations. In both cases, we observe reasonable agreement between data and expectation 
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and the two data sets from the different /3-values scale nicely onto each other, after we have 
corrected for tree-level lattice artefacts. 

In Fig. |1^, the spin-spin potential V4 is displayed in lattice units for the two f3- values. An 
oscillatory behavior is observed which is similar to that of the lattice 5-function, expected at 
tree-level, Eq. (fn|). Moreover, the two data sets nearly coincide with each other, in distinct 
violation of scaling. Corrections to the 5-function, which might scale with an appropriate 
dimension, should account for the differences between the two data sets at small R. 



D. Momentum-dependent potentials 



We intend to compare the MD potentials to Eqs. ( j36f )-fl87D . Since, in accord with these 
expectations, the MD potentials are rather small, compared to the SD potentials, the data 
suffers more from statistical noise, and we do not attempt to perform fully independent fits. 
In addition, we neglect running coupling effects that have been parameterized by / in the 
case of Vo, and V3. We have to subtract the self-energy related constants C b and C d from 
the data points on V b and Vd, respectively, prior to scaling the data sets onto each other. 
We also correct V b and V c for tree-level lattice artefacts. For this purpose we define, 

Ait ~ 2 
6V b (R) = — V6 itree (R) - — , (128) 

5V C (R) = J^V Citree (R) + ^. (129) 



The tree-level lattice expectations for V b and V c can be computed from Eqs. ( [52] ) and 
We fit the data to the following parametrizations, 

V b (R) = || - \KR + C b + g b 5V b (R), (130) 

^(R) = - \KR + g c SV c (&), (131) 

V d (R) = -~KR + C d , (132) 
y 

where we have constrained the parameters e and K to the values, obtained from the fit to 
the static potential. 

The resulting parameter values are listed in Table [VII] . In accord with the BBMP 
relation Eq. ( Pq) , we find — 2C b — ACd = Cy (within errors) for both f3- values. The corrected 
potentials, 

H,cont(#) = H(R) - g b 5V b (R) - C b , (133) 
V c ,cont(R) = V C (R) - g c SV c (R), (134) 



as well as Vd — Cd and V e are displayed in Figs. p!4}-[r7|. The expectations Eqs. (p6|) -(p7|) are 
included as well (solid curves). The x 2 /^df values of the above fits are larger than 1 for 
V b and V c , which means that the correction for lattice artefacts of these potentials is not as 
successful as it has been in the case of and V3. This can be understood from the fact 
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that the MD potentials are more strongly affected by the infrared behavior of the gluon 
propagator, such that higher order corrections might be important. V e shows substantial 
lattice artefacts too (Fig. [T?D . In the case of V& the small- R data lies below the curve, 
indicating that the 1/r coefficient 2e/3 might have been overestimated. This effect cannot 
be understood in terms of the tiny f/r 2 correction that has been omitted. However, by 
allowing for a 1/r-term with a coefficient of about e/8 in Vd, the expectation can be brought 
into agreement with the data. Within — V c the 1/r coefficient e/2 appears to be slightly 
underestimated. 

We conclude that the data is in qualitative agreement with the expectations Eqs. Q8T5D- 
fl57|), although a quantitative comparison fails as there are indirect indications that Vd and 
V e might contain small Coulomb-like contributions, in addition to the linear term. 



E. Comparison with perturbation theory 

In Figs. |i~8"H2l|, we focus on the small- R behavior of the SD potentials V 2 ', . . . , V4 and 
the MD potential V c . We show only the (3 = 6.2 results, which are in qualitative agreement 
with those obtained at (3 — 6.0. Besides the data points, the figures include the tree- 
level perturbative expressions of Eqs. (55)— (pip and Eq. fl5"3"D. The normalization constants 



c = Cpds have been obtained from fits to the first seven data points. V^ and V3 are well 
described by these one-parameter fits and deviations of the data from a continuous curve 
can be understood in terms of this lattice expectation. For V4 as well as V c , agreement is 
only achieved on a qualitative level. The fit parameters are displayed in Table |VII]| . 

From the analysis of the static potential, we expect c = e — h ~ 0.25, compared to the 
tree-level lattice expectations c = 0.106 and c = 0.102 for (3 = 6.0 and (3 = 6.2, respectively, 
determined from the lattice coupling a s = 3/(2tt/3). In agreement with the perturbative 
expectation, all fitted q decrease with increasing (3. We find c c to be about 5 times as 
large as the naive tree-level value; this factor reduces to 2.4 in the case of V 2 ' and 1.9 for 
V3 and V4 as the relevant gluon momenta within these potentials are larger and thus more 
perturbative. 

In order to investigate if remaining differences between data points and renormalized 
tree-level expectations can be explained in terms of higher order perturbative corrections, 
we attempt to model running coupling effect. The only additional diagrams that contribute 
to Vo at 0{g i ) on the lattice (and in the continuum) are one-loop corrections to the gluon 
self-energy. The renormalization of the coupling, emanating from these diagrams, has been 
computed on the lattice for on-axis separations of the sources J53],|5|]. One can account for 
this correction by building in a running coupling constant a(q) into the gluon propagator 
of Eq. (0), in momentum space. Instead of attempting to compute the correct lattice sum, 
we model this effect by the corresponding continuum expression, 



a(t) 



1 ' < " - 1 



4nb t 




(135) 



with t = ln(<f/A 2 ), b = bx/b 2 , b = 11/(16tt 2 ), b x = 102/(16tt 2 ) 2 , where we replace q 2 by 
its lattice counterpart q 2 = 4 J2i sin 2 (gj/2) |i8|j55|| . A is a QCD scale parameter that can be 
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related to the usual schemes via perturbation theory [|4],^(|^7|]. The difference between the 



correct one- loop lattice expression of Ref. |[54| and Eq. ( 135 ) with b = is small. 

In the continuum, contributions that appear in addition to a pure renormalization of the 
gluon propagator can be resummed into a single running coupling, using renormalization 
group arguments. On the lattice rotational invariance is broken and the direction of q enters 
in addition to its absolute value; hence such arguments do not apply. Bearing this in mind, 
we will nonetheless attempt to model higher order perturbative effects by the continuum 
running coupling of Eq. ( |135|) . 

In the case of the SD potentials V 2 ', . . . , V4 not only the gluon self-energy contributes 
to 0(g A ) but also exchange diagrams between the ears, incorporating a three gluon vertex. 
In the continuum, these can be resummed into an effective running coupling. Due to this 
resummation, the scale parameters Aj (for V i ) can differ from each other. However, one 

finds a 2 = A y mm. 

To remove the unphysical pole at q = A, an infrared protection can be built into the 
propagator by substituting t by = ln(g 2 /A 2 + d 2 ) with a constant d. The smallest mo- 
mentum on the lattice is q = n/(aL a ). We choose d 2 = max(0, e — n 2 /(aL a A) 2 ), where e is 
the Euler constant, to guarantee t > 1; d 2 is negligible at large momenta q ~ 1/a. Notice, 
that within the SD potentials the infrared region is suppressed by powers of q, such that 
the results are rather robust with respect to the choice of d or other specific details of the 
protection scheme. 

Fits of the one- and two-loop running coupling improved expressions to the first 4-8 
data points of each potential have been performed. A is the only free parameter within 
these fits. The results of the one-loop fits to 7 data points are included into Figs. |I~8|-|2"0| (full 
circles). The A-parameters remain stable against the variation of the fit range within errors. 
Since the data is described by the tree-level formulae equally well, we are unable to decide 
at present whether the deviations between expectation and data for V4 can be explained 
entirely in terms of such higher order perturbative corrections. 

In Tables |X| and [X], results on one- and two- loop estimates of A-parameters are presented. 
We observe scaling between the two sets of A-parameters obtained at /3 = 6.0 and j3 = 6.2. 
The two-loop values are about twice as large as the corresponding one-loop values. However, 
the (one- and two-loop) aj(g)-values at a scale q = ir/a are consistent with each other. We 
conclude that the one-loop A-values should be considered as effective and not physically 
meaningful. 

From our one-loop fits to V%, we obtain ay(7r/a 6 . ) = 0. 13 1^4 and ay(7r/a 6 2 ) = 0.124l| at 
(3 = 6.0 and (3 = 6.2, respectively. The corresponding two-loop results, ay(7r/a 60 ) = 0.128tg 
and oiy(i{ / — 0. 121^5, are in nice agreement with these numbers, while from the average 
plaquette |47| we obtain ay(7r/a 6 . ) = 0.149 and ay(7r/a 6 . 2 ) = 0.138. We conclude that 
higher order perturbative corrections as well as 0(a 2 ) discretization errors, under which 
localized quantities like the plaquette are more likely to suffer, are responsible for the 2<j-3cr 
deviations between the ay-values, determined from two different observables. 
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VI. APPLICATION TO QUARKONIA SPECTRA 



With the potentials derived from quenched QCD we would now like to predict exper- 
imental quarkonia levels. The spectrum will be computed numerically from the two-body 
Hamiltonian, the structure of which will be summarized in the next subsection. 



A. The Hamiltonian 

Within the spectroscopy study, we restrict ourselves to the equal mass case m = mi 
7T7.2- The starting point is the Hamiltonian, 



H 



5m) +H + 5H kin + 5H CC + 5H sd + 5H md , 



(136) 



where 



m 



+ 



V(r)-5H cc (r) 



(137) 



contains the Coulomb-like part within the relativistic correction to the central potential. We 
numerically solve the radial Schrodinger equation for Ho, and treat 



5 H } 



kin 



5H 

4m 3 ' 



c 3 (m) 



ire5 3 (r) 



5H 



in 



sd 



sd j 



5H 



md 



V, 



md 



(138) 



as perturbations pT |. 

For the particular parametrizations Eqs. (|8~2"D and (|33"D, one obtains the central potential 
[Eq. (0)], 



V(r) 



KT 



(Wni) c 3 (m 2 J 



Till 



m 2 



(2k -b)- +ATie5 3 (r) 
r 



(139) 



The perturbation SH CC is due to the last term within this equation. We have omitted the 
constants Cy, and from the above formula. The latter two of these contributions 
cancel each other within the statistical accuracy of our lattice results while, as we shall see 
below, Cy can be absorbed into a redefinition of the quark masses. The n/r- and <5 3 (r)-terms 
have their origin in the Darwin interaction while the 6/r-term is due to V 2 V^ E . The mass 
dependent correction terms explain the phenomenological flavor dependence of the central 
quarkonium potential, as obtained from fits to the spin-averaged charmonia and bottomonia 
spectra J58],(§. 

Two-particle bound states can be classified by a radial excitation n, the orbital angular 
momentum L, the total spin S = 0, 1 (S = Si + S 2 ), and a total angular momentum 
J = L — S, L, L + S (J = L + S). Conventionally, the states are labeled by n 2S+1 Lj 
where the letters S, P, D, F are used for L = 0,1, 2, 3, respectively. From parametrizations 
Eqs. (|82|), (184]) and (|35l), we find (for equal masses), 
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V sd (r) 



' _K Ac 2 {m){e-h) - e \ Lg 
r 2r 3 / 



(140) 



+ 7 (7t|(m) - 3) 87r(e - /i)5' 



'r 



with 



S1S2 

3 

LS 



2 



L(L + 1)- 5(5 + 1)), 
6(LS) 2 + 3LS - 25(5 + 1)L(L + 1) 



6(2L - l)(2L + 3) 



(141) 
(142) 
(143) 



The one-loop values of the coefficients 02(777) and 03(771) for m = and m = m c at our 
lattice spacings can be found in Table |TJ. The values of the parameters e, h and b are listed 
in Tables [V] and [V]. 

])-(H3), we find the MD correction [Eq. (0)], 



Based on the parametrizations Eqs. 

JC S 



where 



Knd(r) 



K 



6r 



/C 
"6 



_£_ 

2^2 



zrp, 



1 1 

+ 



2 



ml 



mirri2 



£ 



771x1712 



(144) 



(145) 



i.e. /C is a dimensionless parameter while £ carries the dimension m 2 . Note that a string 



of constant longitudinal electric field with energy density k 5S , connecting two pointlike 



particles with masses mi and m 2 , gives rise to the classical correction term — §rL 2 , which 



6 r' 



appears in the above V m &. One obtains the following expectation value within wave functions 
that obey the Schrodinger equation, 



(146) 



such that Eq. (H3) can readily be treated as a perturbation. 

We have neglected the constants Cd and C& of Vd and V& from Eq. ( |144| ). Inclusion of 
these terms would result in a correction, 



AV r 



md 



2 "I 2 ^ 

mf 777-2 / 



777! 7772/ 



V 



(147) 



where we have exploited the relation Cy = —2Cb — 4Cd from the BBMP constraint, Eq. (p7|). 
Under the assumption that \Cf,\ <C 2|C d | (which is supported by our numerical results, the 
MAL picture and tree- level perturbation theory), the above shift in the MD potential can 
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be absorbed into a redefinition of the quark masses of Eq. fll6|) at the given order of the 
Hamiltonian: -> m ; + Cy/2. Notice, that the parameter m of Eq. ( |136| ) differs from 
that of Eq. ([H]) by this constant. If we interpret the fit parameter Cy, obtained at a lattice 
spacing a, as the self-energy of a static quark, E^(fi = ir/a) ~ a~ 1 Cy/2, the combination 
m — 5m within Eq. ( |136[ ) should approach the heavy quark pole mass. 



B. Spectroscopy results 

The physical scale, the quark mass m and the energy shift 8m have to be fixed from 
experiment, before predictions can be made. Note that the parameter m within the Hamil- 
tonian Eqs. (|136|) -( |i~45| ) is not the bare quark mass, but contains the static quark self-energy. 
The value of the dimensionful parameter k determines the scale. We attempt to estimate 
the (presumably small) impact of the parameter 5m on our results. Therefore, we follow two 
strategies: we minimize the squared differences between our predictions and experimental 
levels under the assumption 5m = 0, with respect to m and k (method 1). Alternatively, we 
determine m and k from minimizing deviations from the splittings M(n +1 Lj) — M(l 3 <Si). 
Subsequently, 5m is tuned to reproduce the 1 3 S\ experimental state (method 2). The latter 
method results in the ratios 5m^,jm^^ m 0.04 and 5m c /m C 2 ~ 0.22. However, for ratios 
of the scales, determined by use of the two methods, we obtain (n 2 / Ki) 1 ^ 2 ~ 1.0025 and 
(^2/ ' Ki) 1 / 2 ~ 1.035 for bottomonia and charmonia states, respectively. The tiny size of the 
deviations of these ratios from unity can be understood from the fact that a simple rescaling 
of the mass works rather well as the following ratios indicate: (m 2 — 5m)jm\ « 1.0015 
and 1.012, again for bottomonia and charmonia, respectively. This illustrates that although 
5m can be quite substantial, and thus the uncertainty in the quark mass can be large, the 
impact of this parameter on the predicted spectrum is negligible. Hence, we adopt method 
1, i.e. we set 5m = and allow for two free parameters, m and k. 

In the case of bottomonia, ten states have been observed and we minimize our data with 
respect to all these states. The BB threshold is at about 10.55 GeV. For charmonia, we 
chose to optimize the spectrum only with respect to the seven states below the DD threshold 



at about 3.7 GeV. The results for (3 = 6.0 and (3 = 6.2 are displayed in Tables |XI] and [XII 
for bottomonia and charmonia, respectively^]. We find agreement on the level of about 2-3 
MeV between the results obtained at these two lattice spacings for bottomonia, compared 
to about 10 MeV for charmonia. The differences are likely to reflect the uncertainties in the 
matching coefficients C2{m) and c^{m) which increase with decreasing quark mass. 



The (3 = 6.2 results are compared to experiment in Figs. 22 and E3l In order to estimate 



the effect of quenching, we have included the results obtained for the parameter value e 



0.40, which is our estimate, based on Ref. [14], of the value one might obtain with three 



active flavors of sea quarks. In the case of bottomonia states, an average deviation between 
prediction and experiment of 12.4 MeV at /3 — 6.0 and 12.8 MeV at f3 — 6.2 is observed. 



The effect of the statistical errors of the fit parameters on the spectrum is negligible, in com- 
parison to the systematic uncertainties of the approximation, particularly those of the matching 
constants Ci(fi,m). For this reason, we do not attempt to include any errors into the tables. 



30 



With a value e = 0.40, this deviation is reduced to 9.5 MeV. For charmonia, we obtain an 
average deviation of 22.0 MeV for both /3-values. The parameter choice e = 0.40 changes 
this to 23.0 MeV, indicating that the charmonium spectrum is rather insensitive towards 
quenching effects on the running of the QCD coupling. This can be understood from the 
fact that the wave functions are broader, such that the spectrum is less affected by short 
distance physics. 

From our fit to the bottomonium spectrum, we obtain the following parameter values, 
both at (3 = 6.0 and (3 = 6.2, 

= 468 MeV, m b = 4.68 GeV. (148) 

The above string tension value yields a Sommer scale 0] tq ~ 0.49 fm. tq denotes the 
distance at which the condition r 2 dVo/dr = 1.65 is satisfied. With e = 0.40, we find 
y/7i = 452 MeV and mb = 4.72 GeV. The scale ro remains unaffected under this change in 
e. From the charmonium spectrum, we find, 

= 450(4) MeV, m c = 1.33(1) GeV. (149) 

The errors correspond to the variation between the results obtained at the two /3-values. In 
comparison, the sea quark model with e = 0.40 yields ps 440 MeV and m c w 1.38 GeV. 
The above results are consistent with pole masses ml° le = 4.7(2) and m^ ole = 1.4(2) GeV, 
where the errors are estimated from the <5m/m-ratios. 

From the fit to ten bottomonium states, we find lattice spacings a -1 ps 2.1 GeV and 
a -1 ps 2.9 GeV for the two /3-values, respectively, which are in reasonable agreement with 
estimates from the light hadron spectrum. If we forced the average of the 2 3 S — 1 3 S and the 
1 3 P — 1 3 S splittings to coincide with the experimental counterpart, as it is normally done 
in NRQCD studies, we obtain a -1 ~ 2.5 GeV and a -1 ~ 3.4 GeV, respectively, which is in 
agreement with NRQCD estimates M. As a result, however, the 2 3 P masses would come 
out to be significantly heavier than in experiment. 

In order to investigate the reliability of the nonrelativistic approximation, we have com- 
puted average radii and velocities of various quarkonia states. The results are displayed 
in Table |XIII| . The heavy quark velocities within bottomonia range from {v 2 ) = 0.07 to 



(vl) = 0.11, while for charmonia we obtain the interval, 0.27 < (v 2 ) < 0.52. The radial 



bottomonia wave functions g n i(r) [ip n i m (r, O) = [g n i( r )/ r ] Yim(ty] are shown in Fig. |4. 

We attempt to estimate the approximate size of 0(v 4 ) corrections, using the ratio R = 
(v 4 )/(v 2 ) and find, Rb ~ 0.1 and R c ~ 0.4. Under the assumption that the coefficients 
of such corrections have the same size as those of the 0(v 2 ) corrections, we estimate an 
uncertainty of 4 and 15 MeV for bottomonia and charmonia, respectively, due to neglecting 
higher order terms in v. The uncertainty of the matching coefficients between QCD and the 
effective nonrelativistic theory is another source of systematic biases. We have assumed that 
the coefficient Ci(m), in front of the correction to the kinetic energy, equals 1. As can be 
seen from Table |TTTJ, such coefficients can easily differ by as much as 20 % from this tree-level 
value in the case of bottomonia and by a factor 2 for charmonia. Such an effect on c\ (m) 
would result in shifts of certain bottomonia states by about 4 MeV and charmonia states of 
up to 50 MeV. Higher order corrections to c<i{jn) and c^{m) will also have an effect but at 
present the value of ci(m) constitutes the dominant uncertainty. 
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We conclude that the deviations between experiment and prediction can be explained 
in terms of quenching and higher order relativistic corrections, i.e. 0{y 4 )-terms as well as 
0(lnm//i) and Oimj /i mm/ /z) uncertainties in the matching coefficients of the 0(v 2 ) terms. 
Inclusion of sea quarks seems to improve agreement with experiment but is unlikely to 
reduce deviations by more than an average of 4-5 MeV per state. In the case of bottomonia, 
we estimate the impact of higher order correction terms to be about twice as large. For 
charmonia the effect of 0(v 4 ) corrections might be as large as 10-20 MeV while the impact 
of the matching constants of the 0(v 2 ) terms is even larger. Thus the agreement on a 20 
MeV level appears to be somewhat fortuitous. However, this outcome is not a complete 
surprise since many effects seem to affect the spectrum as a whole, rather than individual 
splittings, and can compensate each other. 

In this first glimpse at the spectrum, we have not yet included a running coupling into 
the parametrization of the potentials. The SD and MD corrections as well as the correction 
to the kinetic energy have so far been treated as first order perturbations only. We will 
improve on these two points in a detailed spectroscopy study f61fl , in which we are going 
to elaborate on the effect of higher order relativistic uncertainties on individual states in a 
more systematic manner. 



VII. CONCLUSIONS AND OUTLOOK 



We have determined the complete 0(v 2 ) relativistic corrections to the static interquark 
potential in SU(3) gauge theory. We find reliable renormalized potentials with good scaling 
behavior. As in our SU(2) study [^T| we report clear evidence for a 1/r 2 scalar exchange 



contribution in the long range spin-orbit potential V[ at the level of 20 % of the Coulomb 
part of the static potential at inverse lattice spacings of 2-3 GeV. The other SD potentials 
are found to be short ranged and are well understood by means of perturbation theory. 
From V 2 ', we obtain the result «y /= V) = 0.124 ± 0.005 ± 0.003 at \i « 9.2 GeV, where the 
first error is statistical and the second one accounts for the differences between one- and 
two- loop estimates. This value is significantly smaller than the estimate ay(/i) = 0.138, 
obtained from the average plaquette. 

All MD potentials contain contributions, that are linear in the quark separation, and 
are in qualitative agreement with minimal area law expectations. The potential V 2 V^, 
that modifies the central force, is found to be Coulomb-like and has a significant effect on 
spectroscopy since it increases the effective Coulomb force by 2 % in the case of bottomonia 
and by as much as 35-40 % for charmonia. A similar behavior is expected from dual 
QCD@. 

As an application, quarkonia spectra are determined. We are able to reproduce the ex- 
perimental levels with an average error of 12.5 MeV for T states and 22 MeV for J/ip states. 
A reduction of these deviations should be achieved by incorporating improved parametriza- 
tions of the lattice potentials, that account for a weakening of the effective QCD coupling at 
small separations, into the spectroscopy. Such a refined analysis is in progress |nj . We esti- 
mate a further improvement of up to 4 MeV per state if dynamical sea quarks are included, 
while higher order relativistic corrections and uncertainties in the perturbative matching 
constants might shift certain levels by as much as 5-10 MeV for bottomonia- and up to 
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50 MeV for charmonia-states. Uncertainties in the (perturbative) matching constants be- 
tween QCD and the effective nonrelativistic theory are likely to have a larger impact than 
0{v A ) corrections, and should be reduced. Our results are compatible with heavy quark pole 
masses m p b ole = 4.7(2) GeV and mP° le = 1.4(2) GeV. 

The approach presented in this article can be used to obtain optimized wave functions for 
creation of a quarkonium state with particular quantum numbers within the complementary 
lattice NRQCD method |63|| . We intend to extend this application to discrete finite boxes 
with periodic boundary conditions, in order to shape even better basis states and to simulate 
finite size effects that one might expect in lattice NRQCD studies. From Fig. it is obvious 
that on volumes with a spatial extent of typically less than 2 fm, excited state wave functions 
become squeezed, and the corresponding energy eigenvalues might be significantly affected. 

Application of the Schrodinger-Pauli approach to the spectrum of B c states as well as a 
determination of quarkonia decay constants is in progress. It appears worthwhile to consider 
calculations on anisotropic lattices, to reduce systematic uncertainties on the potentials, 
arising from the temporal discretization of the lattice. In order to keep uncertainties in the 
perturbative matching constants between the effective theory at a scale fi = it /a and QCD 
small, one would like to operate at spacings a ~ ir/m where m might be either the bottom 
or the charm quark mass. The latter would require an improved lattice action. 
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TABLES 



TABLE I. Simulation parameters. The physical scale has been obtained from yfR, = 468 MeV. 



(3 = 6.0 [3 = 6.2 

V = L\L T 16 3 32 4 

a/fm 0.092 0.067 

a~ l /GeV 2.14 2.94 

L a /fm 1.47 2.15 

n con f 420 116 



TABLE II. Renormalization constants for magnetic and electric ears, compared to their tadpole 
estimates, Z tadpole = 1/U a - 

J3 Un -^tadpole Ze 

6.0 0.593682(5) 1.6844 1.6777(2) 1.6216(4) 

6.2 0.613631(3) 1.6296 1.6249(1) 1.5782(1) 



TABLE III. Matching constants between QCD and the effective Hamiltonian, Eqs. (|l6|)-(^ll), 
for bottom and charm quark masses at (5 = 6.0 and (3 = 6.2. 



C2 C3_ 

m b , P = 6.0 1.034 1.181 

m b , (3 = 6.2 1.065 1.344 

m c , /3 = 6.0 1.220 2.159 

m c , (3 = 6.2 1.257 2.353 



TABLE IV. Fit parameters to the static potential, Eqs. (p5|)-(p7|). 



Parameter (3 = 6.0 (3 = 6.2 Average value 

K 0.0479(12) 0.02536(35) K 

e 0.324(17) 0.321(7) 0.321(6) 

/ 0.042(12) 0.051(5) 0.0082(8)/VT 

Cy 0.6648(78) 0.6404(26) 

g 0.301(4) 0.252(2) 
X 2 /N DF 0.4 0.6 
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TABLE V. Parameter values from fits to relativistic corrections to the static potential that 



are relevant for spectroscopy [Eqs. (118), ( |119| )1. The error in square brackets is the systematic 
uncertainty. Where not stated separately, it has been included. 



Parameter 



(3 = 6.0 



= 6.2 



Average value 



b 
h 

e — h 



0.150(21)[01] 
0.071(5)[11] 
0.253(6)[11] 



0.097(13) [03] 
0.065(3) [8] 
0.256(7) [8] 



3.36(36)K 
0.067(9) 
0.255(10) 



TABLE VI. Fit parameters for the SD potentials V£ and V 3 [Eqs. (g2§-(|2^)]. 



Parameter 


(3 = 6.0 


(3 = 6.2 




0.274(34)[17] 


0.239(16)[58] 


h 


0.023(32) [07] 


0.061(15)[05] 


92 


0.138(40)[01] 


0.133(19)[06] 


xI/Ndf 


1.0 


1.1 




0.253(15)[33] 


0.230(45) [28] 


h 


0.054(11)[41] 


0.047(36) [36] 


93 


0.165(21)[34] 


0.171(07)[29] 


xI/Ndf 


0.4 


1.0 



TABLE VII. Fit parameters for the MD potentials V b , V c and V d [Eqs. (pC| ), (^32|)]. 



Parameter 


(3 = 6.0 


(3 = 6.2 


c b 


-0.0824(71) [2] 


-0.0681(38) [1] 


9b 


0.196(49) [2] 


0.156(41)[1] 


xI/Ndf 


3.0 


5.9 


9c 


0.304(35) [4] 


0.219(12)[2] 


xI/Ndf 


2.4 


1.9 


c d 


-0.116(4)[37] 


-0.122(2)[32] 


X 2 JNdf 


0.5 


1.2 
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TABLE VIII. The constant c = Cfoi s from the weak coupling analysis. 



Potential 


p 


c 


X 2 /N DF 


V c 

vi 
v 3 
v A 


6.0 
6.2 
6.0 
6.2 
6.0 
6.2 
6.0 
6.2 


0.58(9) 
0.48(4) 
0.251(19) 
0.238(12) 
0.192(22) 
0.182(15) 
0.217(17) 
0.199(13) 


4.1 
5.2 
0.84 
1.30 
0.34 
0.46 
10.3 
12.0 




TABLE IX. A-parameters from the 


one-loop running coupling analysis. 




Potential 


P 


A/v^ 


X 2 /N DF 


vi 
v 3 
v A 


6.0 
6.2 
6.0 
6.2 
6.0 
6.2 


0.1812 
0.20^ 
0.090l|| 
0.09lllg 
0.08711? 
0.0561H 


1.25 

0.93 

0.59 

1.15 

4.5 

12.7 


TABLE X. A-parameters from the two-loop running coupling analysis. 


Potential 






X 2 /N DF 


n 

v 3 
v A 


6.0 
6.2 
6.0 
6.2 
6.0 
6.2 


0.411" 
0.451? 
0.22ljj 
0.22l^ 
0.19±| 
0.1311 


1.32 

0.96 

0.59 

1.15 

5.5 

14.2 
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TABLE XI. The bottomonium spectrum. 





P = 6.0 


P = 6.2 


e = 0.40 


experime 




9.477 


9.476 


9.415 




l 3 S l 


9.521 


9.526 


9.504 


9.460 


2 1 S 


9.980 


9.980 


9.961 




2 3 S X 


10.007 


10.010 


10.008 


10.023 


3^0 


10.328 


10.328 


10.311 






10.351 


10.354 


10.348 


10.355 


4 1 S 


10.619 


10.619 


10.597 




4 3 5i 


10.640 


10.642 


10.630 


10.580 


l'Pi 


9.879 


9.879 


9.889 




1 3 P 


9.866 


9.866 


9.867 


9.860 


l 3 P l 


9.878 


9.878 


9.886 


9.892 


1 3 P 2 


9.882 


9.883 


9.895 


9.913 


2 1 P 1 


10.238 


10.238 


10.243 




2 3 P 


10.226 


10.225 


10.223 


10.232 


2 3 Pi 


10.237 


10.237 


10.240 


10.255 


2 3 P 2 


10.241 


10.242 


10.249 


10.269 


l l D 2 


10.120 


10.121 


10.136 




l 3 D 1 


10.121 


10.121 


10.134 




l 3 D 2 


10.121 


10.122 


10.137 




1^3 


10.119 


10.120 


10.137 
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TABLE XII. The charmonium spectrum. Only the states that are identified by an asterisk lie 
below the DD threshold and have been used to fix the scale k and quark mass m c . 



(2S+l) Lj 


= 6.0 


P = 6.2 


e = 0.40 


experiment 


1% 


3.010 


3.001 


2.966 


2.980* 


1 3 S 1 


3.134 


3.143 


3.175 


3.097* 


2 1 S 


3.591 


3.582 


3.560 


3.594* 


2 3 5i 


3.685 


3.688 


3.705 


3.686* 


3% 


4.017 


4.004 


3.978 




3 3 Si 


4.102 


4.098 


4.106 


4.040 


4^0 


4.371 


4.354 


4.324 




4 3 5i 


4.452 


4.442 


4.442 


4.415 


I'Pi 


3.468 


3.472 


3.486 




l 3 Po 


3.452 


3.451 


3.442 


3.415* 


l 3 Pi 


3.479 


3.482 


3.490 


3.511* 


1 3 P 2 


3.465 


3.471 


3.491 


3.556* 


2 1 P 1 


3.915 


3.913 


3.916 




2 3 P 


3.893 


3.886 


3.870 




2 3 Pi 


3.922 


3.919 


3.917 




2 3 P 2 


3.916 


3.914 


3.924 




l l D 2 


3.764 


3.767 


3.782 


3.770 


1 3 P>! 


3.791 


3.790 


3.796 




1 3 £> 2 


3.777 


3.778 


3.792 




1 3 P> 3 


3.744 


3.748 


3.770 





TABLE XIII. Average velocities and quark separations in bottomonia and charmonia. The 
values correspond to the parameter choice e = 0.40, which has been used to model the effect of 
dynamical sea quarks. 



15 0.080 0.27 0.24 0.43 

2S 0.081 0.35 0.51 0.85 

3S 0.096 0.44 0.73 1.18 

4S 0.112 0.52 0.93 1.47 

IP 0.068 0.29 0.41 0.67 

2P 0.085 0.39 0.65 1.04 

ID 0.075 0.34 0.54 0.87 
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FIGURES 



R 




FIG. 1. Lattice definition of the nominator of Eq. ( [Iq ) for the example of F± being an electric 
ear and F2 being a magnetic ear. 




FIG. 2. Corrected static potential Vb )CO nt at /? = 6.0 and (3 = 6.2. The fit curve corresponds to 



the parametrization Eq. (117) with the parameter values listed in the last column of Table IV 



41 



1.04 



1.035 - 




FIG. 3. The ratio of tadpole and HM renormalization constants Q [Eq. ( |110| )1 as a function of 
t at R = 4, p = 6.2. 



O 




FIG. 4. The ratio Q as a function of R at (3 = 6.2 (for large t). 
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FIG. 5. Test of the Gromes relation, Eq. (|2q) . The combination V% — V{ is compared to the 
static force as obtained from the parametrization Eq. (117). 
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FIG. 6. Test of the BBMP relation, Eq. (|2q). The combination Vj, + 2Vd is compared to its 
expectation as obtained from the parametrization Eq. (|117|) of the static potential. 
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FIG. 7. Test of the BBMP relation, Eq. (]27|). The combination V c + 2V e is compared to its 
expectation as obtained from the parametrization Eq. (|117| ) of the static potential. 
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FIG. 8. The potential V 2 ^, together with a fit curve of the form V 2 V r Q E (r) = — b/r, with 
b = (0.86 ± 0.05 GeV) 2 . The constants have been subtracted from the data points. 
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FIG. 9. The potential V 7 a at (3 = 6.2 in lattice units (statistical errors only). 
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FIG. 10. The spin-orbit potential V{, together with a fit curve of the form —V[(r) = k + h/r 2 . 
with h = 0.067(9). 
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FIG. 11. The spin-orbit potential V^'ccmt ^ comparison to the continuum expectation from 
Eq. Q. 
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FIG. 12. The spin-spin potential V3 :Con t in comparison to the continuum expectation from 
Eq. O. 
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FIG. 13. The spin-spin potential V4 for the two /3-values in lattice units. 
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FIG. 14. The MD potential H.cont in comparison to the continuum expectation from Eq. 
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FIG. 15. The MD potential — V CC ont m comparison to the continuum expectation from Eq. (|86| 
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FIG. 16. The MD potential — in comparison to the continuum expectation from Eq. (87) 
The constants Cd have been subtracted from the data points. 
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FIG. 17. The MD potential —V e in comparison to the continuum expectation from Eq. (| 
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FIG. 18. Comparison of the lattice potential at /3 = 6.2 to tree-level lattice perturbation 
theory [Eq. ([58])] and to the one-loop model of Eq. (135). 
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FIG. 19. Same as Fig. |18| for V 3 [Eqs. Q5g) and flKJ)] 
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FIG. 20. Same as Fig. for V A [Eq. © 
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FIG. 21. Comparison of the lattice potential — y c at /? = 6.2 to tree-level lattice perturbation 
theory, Eq. d65 
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FIG. 22. The bottomonium spectrum. The e = 0.32 results are from the f3 = 6.2 analysis. 
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FIG. 23. The charmonium spectrum. The e = 0.32 results are from the (3 = 6.2 analysis. 
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FIG. 24. Bottomonium wave functions. 
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